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Lattice Boltzmann simulations have been very successful in simulating liquid-gas and other multi- 
phase fluid systems. However, the underlying second order analysis of the equation of motion has 
long been known to be insufRcient to consistently derive the fourth order terms that are necessary 
to represent an extended interface. These same terms are also responsible for thermodynamic 
consistency, i.e. to obtain a true equilibrium solution with both a constant chemical potential and a 
constant pressure. In this article we present an equilibrium analysis of non-ideal lattice Boltzmann 
methods of sufficient order to identify those higher order terms that lead to a lack of thermodynamic 
consistency. We then introduce a thermodynamically consistent forcing method. 



I. INTRODUCTION 

The standard lattice Boltzmann approach leads to an 
ideal gas equation of state. Several different approaches 
to simulate non-ideal fluids with lattice Boltzmann have 
been introduced in the past. The main application has 
been to simulate phase separation, although other appli- 
cation like an increase in the speed of sound have also 
been considered. There are two different philosophies to 
introduce the non-ideal terms. The first is guided by an 
atomistic picture and local interactions are introduced 
through a forcing term P, |3| . The second starts from the 
Navier-Stokes equation for a non-ideal gas and tries to 
match the hydrodynamic limit of the lattice Boltzmann 
equation with this macroscopic equation 0, Q . 

Whatever the underlying philosophy, each approach 
leads to some non-ideal equation of state. Knowing 
the equation of state is sufficient to predict the phase- 
behavior from equilibrium thermodynamic arguments. 
The equilibrium densities are determined through the 
Maxwell construction Many approaches fail this test 
and the resulting phase-diagrams deviate from the theo- 
retical one. 

For many simple applications where phase-separation 
is only required to form a nearly immiscible system, these 
thermodynamic details may be of limited importance. 
For simulations of phase-separation, however, such de- 
tails can be crucial. In this paper we will examine why 
lattice Boltzmann approaches using a force to introduce 
the non-ideal equation of state fail to obtain the cor- 
rect thermodynamic behavior. We explicitly identify the 
terms that lead to the non-thermodynamic behavior of 
these methods. 

The paper is organized as follows: we first identify the 
hydrodynamic equations that we intend to simulate. We 
then discuss the equilibrium behavior of these equations. 
We then introduce a general lattice Boltzmann method 
which introduces non-ideal terms through either forcing 
or pressure terms. The hydrodynamic limit of this ap- 
proach is presented. This allows us to identify how we 
can incorporate the non-ideal pressure contributions by 
either incorporating a bulk force or by altering the pres- 
sure moment of the equilibrium distribution. 

These two different methods for non-ideal systems are 



equivalent to the old Y. Chen method or the more 
recent extension by He et al. |^ 01 when we use the bulk 
force and the Holdych correction to the Swift model Q 
when we directly alter the pressure moment. However, 
we are not able to consistently introduce surface tension 
effects at this point since those appear as higher order 
derivatives that are not derivable by a second order ex- 
pansion 0. 

These higher order derivatives, however, are necessary 
to achieve thermodynamic consistency. We introduce a 
near equilibrium analysis of sufficient order to consis- 
tently derive these higher order gradient terms. This 
analysis uncovers correction terms for pressure and forc- 
ing methods. Despite the fact that we only perform a 
5th order analysis we are able to identify the correction 
terms exactly when the fiuid is not advected with respect 
to the lattice. With this knowledge we are then able to 
formulate a forcing method that achieves thermodynamic 
consistency. 

II. EQUATIONS OF MOTION FOR A 
NON-IDEAL GAS 

As a simple example for a non-ideal fluid we will exam- 
ine a single component fluid that can undergo a liquid-gas 
phase-transition. For simplicity an isothermal system is 
considered. Such a system has an underlying free energy 
of the form 

J[f{p) + IiVp,\7^p,...)]dx (1) 

where p is the density, f{p) is the bulk free energy and 
I{Vp, V^p, . . .) is a gradient expansion of the interfacial 
free energy. The lowest order term for the free energy 
is ■|(V/9)^ and usually only this term is considered. The 
equations of motion for a non-ideal gas are given by the 
continuity equation 

dtp + V(pu) = 0, (2) 

u is the velocity of the fluid, and the Navier-Stokes equa- 
tion of a non-ideal gas 

pdtu + pu.Vu pV/i + Vcr. (3) 
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Here ^ = ^ is the chemical potential of the non-ideal 

gas and a ~ Vu + (Vu)^ — |V.ul is the Newtonian 
stress tensor. 

Let us consider how this system will approach equilib- 
rium. The equilibrium will be time independent and the 
flow will be uniform, i.e. u =const. So the equation Q 
reduces to 

V/i = (4) 

so that the chemical potential will be constant in equi- 
librium. 

Condition Q only guarantees that the chemical po- 
tential is constant in both phases. Bulk thermodynam- 
ics, however, requires that the pressure is also constant 
in equilibrium. This is guaranteed through the thermo- 
dynamic relation 

pV/i = VP (5) 

where P is the pressure tensor. This pressure tensor 
is defined through the properties that it is equal to 
the bulk pressure in the absence of density gradients 
P = pi = [pdpf{p) — f (p)]! and through condition (0) 
in the interfacial areas [l^- With this relation a uni- 
form chemical potential is equivalent to a divergence-free 
pressure tensor. 

As an aside it is interesting to note that this statement 
is more general than bulk thermodynamics. It can even 
be applied to finite systems with curved interfaces. An 
example is an equilibrium drop that will have a constant 
chemical potential, but the pressure inside the drop will 
differ from the pressure outside by the Laplace pressure 
Ap — J, where a is the surface tension. Despite this 
difference in the pressures inside and outside the drop the 
divergence of the pressure tensor vanishes everywhere. 
For such a system with a curved interface the equilibrium 
liquid and gas densities will differ slightly from the bulk 
thermodynamic values . 

The connection between bulk thermodynamics and the 
equilibrium condition pV/i = VP = can be established 
by considering a flat interface between the liquid and gas 
phases. Assume that this interface is orthogonal to the 
X-axis. In this case VP = dxPxx^x is just the derivative 
of a scalar function. Therefore it follows from VP — 
that in the bulk phase pi = Pg- Therefore the solution 
of the differential equation VP — pV/i = for a flat 
interface implies the standard bulk thermodynamic equi- 
librium conditions p; = pg and pi = Pg- 

Let us consider the condition pVp = in some 
more detail. For concreteness' sake let us assume that 
/(Vp, ...) = f(Vp)2. We then obtain p = 9p/(p)-KV2p 
and P = [p9p/(p) - /(p) - npV'p - («/2)(Vp)2]l + 
kV p\7 p. Then the differential equation for a single flat 
interface becomes 

«a|p = dxdpfip) (6) 

subject to the boundary conditions that 

lim dxP^O. (7) 

X — ^±oo 



Because (jBJ is equivalent to both Vp = and VP = 
the limiting values for the density will be the equilib- 
rium densities pg and p;. This conclusion, however, is 
independent of the form of the interfacial energy term 
I{V p, V^p, . . .) and as long as the new differential equa- 
tion has a solution that fulfills the boundary condition 
(|7|) the limiting densities pg and pi will be the same. 

The important corollary of this argument is that while 
there are many possible forms of the chemical potential 
and corresponding pressure tensor that lead to the cor- 
rect bulk phase behavior, arbitrary derivative terms (that 
might arise because of unintentional higher order correc- 
tions to the numerical method) will in general not be 
derivable from an interfacial energy term /(Vp, . . .). In 
such a situation the differential equation can lead to 
bulk densities that do not correspond to the equilibrium 
densities. These situations are the main concern of this 
paper. 

III. LATTICE BOLTZMANN 

The lattice Boltzmann method can be viewed as a dis- 
cretization of the Boltzmann equation. And in the same 
way that the Boltzmann equation describes a gas that 
at long wavelength obeys the hydrodynamic equations, 
the same is true for the lattice Boltzmann method. In 
the lattice Boltzmann method both the space and veloc- 
ity space is discretized and the basic variables are the 
densities fi{x,t) associated with the velocity Vi. The hy- 
drodynamic variables are then the local density 

p(a;,i) = ^/,(a;,i) (8) 

i 

and the momentum 

p{x, t)u{x, t) Mx, t)v,. (9) 

i 

Most lattice Boltzmann methods do not conserve energy 
and instead enforce a constant temperature. One caveat 
is that u is not necessarily the local velocity, as we will 
see below. 

The evolution equation for a non-ideal fluid can be 
written as 

/,(x + Vi,i+l) = /,(x,i)+P,(x,t) + 

i(/0(p) + A,(x,t)-/,(x,t))(10) 

Here is the equilibrium distribution for the ideal gas. 
The Ai represent non-ideal contributions to the pressure 
tensor j^], and Fi are contributions of an external force. 
The force can also be used to mimic interactions in a 
mean-field manner. 

In lattice Boltzmann the Navier-Stokes equations are 
not discretized directly; instead the moments of the equi- 
librium distribution as well as Ai and Fi are chosen such 
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that the momentum moment of a second order expan- 
sion of the lattice Boltzmann equation will give the de- 
sired Navier-Stokes equation Q . A dilemma occurs when 
the pressure tensor itself contains second order derivative 
terms since these terms are formally higher than second 
order in the Navier-Stokes equations. These terms are 
not consistently derived in standard lattice Boltzmann 
expansion techniques. 

So the question arises: how can these terms be con- 
sistently incorporated into a lattice Boltzmann method? 
One may consider simply expanding the lattice Boltz- 
mann equation to higher order, but this will lead to 
Burnett level equations, which is not the desired result. 
The reason that these higher order density derivatives 
appear in the Navier Stokes equation is because the den- 
sity derivatives are not small near an interface and these 
terms are responsible for the surface tension. 

It is difficult and rather cumbersome to extend the 
expansion of the lattice Boltzmann equation to higher 
orders in a general way. So instead we will consider an 
equilibrium (or at least stationary state) interface instead 
and develop a higher order analysis for this situation. 
Doing so uncovers the additional terms that lead to a 
lack of thermodynamic consistency for forcing methods. 
Incorporating these terms allows for thermodynamic con- 
sistency as will be discussed below. 

There are two expansion methods that are regularly 
used to derive the hydrodynamic limit of the lattice 
Boltzmann equations: they are referred to as the Tay- 
lor expansion and the Chapman-Enskog methods. Up to 
second order the methods give identical results, but there 
is debate about the equivalence of the two methods for 
higher order results. In this paper we will utilize the first 
method. 

To establish a starting point for a higher order expan- 
sion we will first present the standard second order ex- 
pansion. We will then show how to choose the moments 
of , Ai, and Fi to simulate a van der Waals gas. 



tribution function. These are given by 

E/" - p (11) 

i 

E/f(v.-u) = (12) 

i 

E/°(v.-u)(v,-u) = pel (13) 

i 

i 

Galilean invariance would require the third order tensor 
Q to vanish. In most standard lattice Boltzmann meth- 
ods this term does not vanish, however, because these 
models contain too small a velocity set. For these mod- 



els vt 



which precludes the presence of the third 



power of u in the ffv'l^ moment. It also fixes the 
temperature to be 6* = I. Therefore most models have a 
Q term given by the third order tensor Q = puuu. The 
effects of the Galilean invariance violation caused by this 
term become noticeable only at large velocities u [ij. 

The non-ideal contributions from the Ai need to con- 
serve mass and momentum and the moments are given 
by 



i 

Y.A^^{v,-n) 

i 

i 

^A^ivia ~ Ua){Vii3 - U/3)(Wi7 - %) 



(15) 

(16) 
Ac^p (17) 

Aaf3Uy + Aa^fUp 



+Ap^Ua. (18) 

The forcing term Fi needs to conserve mass, therefore 



E^.^o. 



(19) 



It also needs to change the momentum by an amount F, 
therefore we choose 



'F 



(20) 



IV. SECOND ORDER EXPANSION 



To obtain the hydrodynamic equations that govern the 
evolution of the slow dynamics of the conserved quanti- 
ties we use a Taylor expansion of (|10|l to second order. 
As usual we will only conserve mass, defined through the 
density © and momentum defined through 

For this iso-thermal model, which does not include a 
conserved energy moment, we require the knowledge of 
the first three velocity moments of the equilibrium dis- 



The second moment of the forcing term is usually taken 
to be zero, but we leave a general term ^E* which we will 
use later to contain a correction term. 



(Vi - u) = 



(21) 



At this point it is worth while to note that the distinction 
between the A and the ^E" terms is somewhat artificial. 
Both terms enter the lattice Boltzmann equation in the 
same way, so that instead of we can insert the term 
A = T\l/ and vice versa. We distinguish between these 
terms to connect our analysis to established methods. 
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Pressure methods l3 only use A and forcing methods 
[1I1I1I1II3 only F and ^. 

Now we need to derive the hydrodynamic limit of the 
lattice Boltzmann equation (|10|l to link this method to 
the hydrodynamic equations (0), lO which we want to 
simulate. 



A. The hydrodynamic limit 

While the Taylor expansion method is in principle well 
known we will present it here again because the higher 
order analysis presented later in this paper uses the re- 
sults and techniques of this approach. The main premise 
of the Taylor expansion approach is that the distribution 
function can be expressed by a Taylor expansion 



f^{x + V^At,t + At) = J2 



(Ml 

k\ 



D''Mx,t). (22) 



where we have defined the derivative operator Di — [dt + 
ViV). For this expansion to be useful we must make 
the assumption that derivatives are small. One could 
phrase the same argument in terms of At, but it will be 
convenient to set At = 1 in (|10|1 . To order our terms we 
will therefore use the derivatives as a small parameter e 
so that 0(9") = 0(e"). We obtain for ^ to second 
order 



/.), (23) 



where we have introduced ff — F^ + Ai. This relates 
non-local fi to local and their derivatives. However, 
the fi remain unknown and we only know the ff in terms 
of the macroscopic quantities. We can use (|23|l to express 
the fi in terms of the equilibrium distribution and higher 
order derivatives: 

/. = /° - tF, - tDJ, + 0{d^) 

= /° - tF, - TD,{f° - tF,) + 0(8'). (24) 

Using this approximation we can now express the lattice 
Boltzmann equation purely as a differential equation in 
terms of the equilibrium distribution and the collision 
term: 

F, + A(/° - tF.) - - - rF,) 

= ^(/°-/0 + O(93)- (25) 

Taking the zeroth order velocity moment - H25|l we ob- 
tain (borrowing the Euler level terms, i.e. the terms of 
0{d), of the momentum equation (|29|l ) the continuity 
equation: 



This leads us to identify the mean fluid velocity as 



1 



u = u F. 

2p 



(27) 



We then obtain the continuity equation 

atP + V(pu) = 0(93). (28) 

Taking the first order velocity moment Vi ()25|l we ob- 
tain the Navier-Stokes level equation: 

pdtU+pu.Vii = -V{pe+A)+F+V(T+VR+0{d^) (29) 

where the Newtonian stress tensor a is given by 

a = iyp{\7u + (Vuf) (30) 

and unphysical terms have been collected in the remain- 
der tensor 

i? = T* - 3v[uy.A + (uy.Af + u.VAl + VQ] + 0(9^). 

(31) 

The kinematic viscosity is given by = (r — ■| )6>. Note 
that while we wrote u in equations (|30|l and Ij^l^ it cannot 
be distinguished from u here because the forcing term is 
0{d) and therefore u — u+0{d). We also note that most 
of the unphysical terms in H31|l violate Galilean invariance 

B. Forcing and Pressure methods for the Non-ideal 



If we set both A and F to zero we obtain the Navier- 
Stokes equation for an ideal gas. To obtain the Navier 
Stokes equation for a non-ideal gas previous research has 
identified two different strategies. One option is to use 
the forcing term to introduce the non-ideal contribution 
to the equation of state 3 . The form presented here was 
first presented in the elegant work of He et al. 0, • In 
this case, which we will refer to as the forcing method, 
we define 



F = - 
* = 0, 
A = 0. 



V.P" 



(32) 



(26) 



where P™'^ = P—nOl is the non-ideal part of the pressure 
tensor. The second approach is based on the idea of 
including the non-ideal pressure in the second moment 
of the equilibrium distribution^,]^. To do this we define 

F = 0, 

^' = 0, (33) 
A = P - pe\ + v{n\/p+ {uVpf + Vi.Vpei), 

where the v terms have been introduced by Holdych et 
al. and later by Inamuro et al. To do this 
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FIG. 1: Comparison of the profile (a) and tlie pliase diagram 
(b) for the van der Waals pressure of equation 13611 obtained 
by the forcing method of equation H82|l with the pressure 
method from equation I33II . There is no unique solution to the 
pressure method. The unique solution of the forcing method 
is not in agreement with the theoretical phase diagram. The 
critical temperature used in (a-c) was 9c — 0.37 at a tempera- 
ture of S = 1/3. In (d) no liquid-gas phase boundaries can lie 
outside the white area and still have equal pressures for the 
liquid and gas phase. 



we have to use the near-equihbrium approximation of 
A = —pO + 0(e), as discussed by Kalarakis [l^ . 

Up to second order in the derivatives both of these 
approaches lead to the same hydrodynamic equations 



dtp + V(pu) 
pdtVL + pu.Vu 







-VP + Vcr 



(34) 
(35) 



For the van der Waals pressure for a critical density of 
1 and a temperature 6 — \ the pressure tensor is given 

by 



P 



3-p 



(36) 



Using this pressure for the two methods we observe 
phase-separation below the critical temperature as shown 
in Figure n There are, however, a number of peculiari- 
ties. 

Let us first discuss the results for the pressure method. 
There is no unique solution. Three different profiles 
for the pressure method are shown in Figure ^ a). All 
of them are stable solutions which are stable against 
small perturbations. Also the density profiles are sharp, 
switching from the liquid to the gas density in only one 
lattice spacing. These profiles correspond to a constant 
pressure profile (Fig. H^b)), which is exact to machine 
accuracy. The chemical potential profiles, however, are 
vastly different as shown in Figure ^c). Only the mid- 
dle profile corresponds to a constant chemical potential, 
and therefore to the thermodynamic equilibrium. In Fig- 
ure^d) we show the phase-diagrams. Only in the white 
area in this diagram can liquid and gas phases have equal 
pressures. We show the results of simulations which are 
initialized with a near constant pressure in the gas and 
liquid phase. The A triangles correspond to simulations 
initialized with the lowest possible constant pressure and 
the V triangles correspond to simulations initialized with 
the highest possible constant pressure. All these simula- 
tions are stable to small perturbations. It appears that 
all states with equal pressures, including but not limited 
to the true equilibrium state, are stable solutions. 

The forcing method leads to a unique profile which 
has an interface extending over several lattice spacings. 
The bulk pressures of the liquid and gas phases agree, 
but there are variations of the pressure H36II around the 
interface in Figure ^^b). The bulk values of the chemi- 
cal potential, however, are not equal. It is therefore not 
surprising that the phase-diagram differs from the true 
equilibrium profile as shown in Figure ^d). 

To understand these results let us recall that expres- 
sion for the pressure (|36|l does not contain any gradient 
terms. This corresponds to /(Vp, V^p, . . .) = in the ex- 
pression for the free energy . This means that we have 
not input any interfacial free energy contributions and we 
therefore expect a sharp interface. But without interfa- 
cial contributions the differential equation VP = is no 
longer a differential equation. Instead this only requires 
that the pressure for both phases is the same, pi = Pg. 
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Any such pressures will fulfill VP = 0, as shown in Figure 
^p. In this Figure numerical solutions at the maximum 
and minimum of the coexistence liquid and gas pressures 
as well as the true equilibrium pressure are shown. Only 
for the true equilibrium densities will /i/ = fig be true 
as can be seen in Figure So clearly VP = pV/i is 
no longer generally valid. There is no guarantee that 
a system without interfacial energy obeying the dynamic 
equations |(5J and (PJ will move towards true equilibrium. 
This explains why the pressure method, which leads to 
the expected sharp interface, can fail to recover the equi- 
librium densities. 

This, however, is not enough to understand the stabil- 
ity of the interfaces to small perturbations. To change 
the liquid and gas densities it is in general necessary to 
move the interface. Because the method leads to a sharp 
interface, moving this interface will lead to states that 
have one lattice point with a density between the gas 
and liquid densities. If the gas density is slightly in- 
creased, the pressure will also increase, favoring a return 
to the original density. Likewise if at one point the liq- 
uid density is reduced this will lead to a lower pressure 
inducing the return to the original density commensu- 
rate with the surrounding pressure. So the stability of 
these non-equilibrium structures occurs because moving 
an interface requires different discretizations of the in- 
terface. And these discretizations lead to position de- 
pendent pressures which counteract the movement of the 
interface. 

The results shown in Figure ^ are at a mean veloc- 
ity of zero. The situation changes when a mean velocity 
is added to the system. In these cases we do not find 
advected stationary solutions, but instead large oscilla- 
tions are observed for sharp interfaces. This violation of 
Galilean invariance is removed by increasing the width of 
the interface, as shown in [T^ . 

In the case of the forcing method the situation is more 
complicated. The forcing method leads to a unique ex- 
tended interface as shown in Figure^. Since we do not 
obtain a sharp interface the pressure P of equation H36() 
is not constant, as seen in Figure^. Therefore higher 
order gradient terms must be present in the pressure for 
this lattice Boltzmann method. These gradient terms 
arise because of higher order terms in the lattice Boltz- 
mann method that were not picked up by our second 
order expansion. The bulk values of the pressure in Fig- 
ure are the same for the liquid and gas phases as is 
to be expected from the condition of mechanical equilib- 
rium. The same, however, is not true for the chemical 
potential in Figure Because the bulk values of the 
chemical potential are not the same the densities also do 
not correspond to their equilibrium values. 

For the pressure method there are two potential reme- 
dies for us to recover the equilibrium bulk densities for 
the liquid and gas phases. The most obvious one, in light 
of the present discussion, is to explicitly include the cor- 
rect gradient terms in the pressure and we will follow that 
route below. A second potential remedy lies in including 
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FIG. 2: The phase behavior for a liquid gas system when the 
explicit interfacial energy is included. For this set of simula- 
tions we used k = 0.1 in the expression for the pressure 1371 . 
All other parameters are the same as for Figure Q 



fluctuation terms in the Navicr-Stokes equations. Incor- 
porating equilibrium fluctuations allows the interfaces to 
move and to select the correct equilibrium bulk densi- 
ties for the pressure method. This approach, however, is 
outside the scope of the current paper. 

While derivatives in the pressure tensor are not consis- 
tent with the second order expansion, it appears reason- 
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able to include the full derivative terms in the pressure 
tensor. This is indeed what was done in the original 
pressure method and since the numerical simulations 
do not indicate the presence of any spurious interfacial 
terms it is not surprising that this approach is successful 
for the Pressure algorithm. 

The usefulness of including the full gradient terms in 
the pressure tensor is less obvious for the forcing method 
where there are clearly already significant, spurious in- 
terfacial terms present that lead to an extended interface. 

So we will now replace (|3ti|l with 



P= p^6»e-K(pV^p+-Vp.Vp) 1 + KVpVp. 

o ~ p o Z 

(37) 

For a numerical implementation the gradient and Laplace 
operator have to be replaced by discrete versions. This is 
problematic because it will in general break the relation 
VP — /oV/i, and thereby also the exact thermodynamic 
consistency. For interfaces that are wide enough though, 
this relation will still hold to good approximation. We 
have not yet been able to identify a discrete derivative 
operator and an expression for the pressure that preserves 
V^P = pV^/z. Such an expression could guarantee that 
a constant pressure is exactly equivalent to a constant 
chemical potential, and therefore it would exactly recover 
equilibrium thermodynamics. 

For the simulations presented in this paper we require 
the one dimensional discrete gradient and Laplace oper- 
ator. We use the discretization 

V^''^f{x)^\{f{x + l)^f{x-l)), (38) 

V2(^V(a:) = fix + 1) - 2/(.t) + f{x - 1). (39) 

For K = the expression leads to But for 

finite k we now expect to find an extended interface for 
the pressure method. We also expect that a constant 
pressure will now reduce the difference in the chemical 
potential for the two phases. The results are shown in 
Figure |21 

For a sufficiently large interfacial energy contribution 
of /v = 0.1 there is now a unique solution for the pressure 
method. This solution also agrees very well with the an- 
alytical phase-diagram. We find that the pressure is con- 
stant up to machine accuracy and the chemical potential 
is very nearly constant. In particular the deviation in the 
bulk value of the chemical potential in the gas and liquid 
phase are less than 3 x 10^"'. 

In the continuous situation the non-definiteness of the 
stationary state is limited strictly to k = 0. For a dis- 
crete approach V^P — pSI^ p. is not generally valid, and 
therefore it is not apriory clear that an arbitrarily small 
interfacial term in the pressure will guarantee a conver- 
gence to the equilibrium solution. As we have seen be- 
fore it is not even guaranteed that it will converge to a 
unique solution. We therefore performed a set of simu- 
lations for different values of n from to 0.1 with ini- 
tial conditions that correspond to the extreme values of 




FIG. 3; Convergence of the pressure in the two-phase system 
for different initial conditions for 6 = 0.37. A unique solution 
is only found for large enough values of k. For a detailed 
discussion of this effect see the main text. 



the equal pressure as well as the equilibrium solution. 
The simulation results are shown in Figure |31 We see 
that non-uniqueness survives for finite k. The range of 
possible pressures is rapidly reduced until a near-unique 
solution is found at k larger than about 0.08. 

The origin of this non-uniqueness even in the presence 
of an interfacial energy lies in the discretization of the 
interface. Moving the interface over the lattice requires 
a re-discretization of the interface. Each of these dis- 
cretizations will have a different interfacial free energy 
X;i(K/2)(V/9)^ To visualize this imagine a sharp inter- 
face of Fig This interface has only two points with 
non- vanishing gradient and Laplace operators. If this in- 
terface moves there will have to be at least one point with 
intermediate value of p and now there are at least three 
points with non-vanishing gradient and Laplace terms. 
Therefore moving the interface will change the total free 
energy. Mass-conservation generally implies that chang- 
ing the liquid and gas densities requires changing the 
relative gas and liquid volumina. Therefore changing the 
densities requires that the interfaces move. But if mov- 
ing the interface requires passing a local maximum in the 
free energy, this implies that the system is at least in a 
metastable state. And a lattice Boltzmann method with- 
out fluctuations cannot escape such a metastable state. 
So to move an interface that is well aligned with the 
lattice one has to overcome a free energy barrier with re- 
spect to the less favorable discretization of the interface 
and the simulation can get stuck in a non-equilibrium 
configuration. 

The Forcing method is also noticeably improved 
through the inclusion of the gradient terms. The differ- 
ence in the bulk chemical potential is about halved from 
0.03 to 0.015. This is still about 50 times larger than the 
difference for the pressure method. This difference in the 
chemical potential translates to a deviation of the liquid 
and gas densities from their equilibrium values. This can 
be clearly seen in Figure In order to understand the 
origin of the large deviations of the forcing method from 
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the thermodynamic equilibrium we need to uncover the 
spurious gradient terms in the effective pressure. To do 
this we need to improve our expansion method to consis- 
tently derive those derivative terms. 



V. HIGHER ORDER NEAR EQUILIBRIUM 
EXPANSION 

To identify the higher order terms of the equilibrium 
structure we perform a 5th order expansion around an 
equilibrium profile which is stationary, but possibly ad- 
vected with a constant velocity U . As before we make 
the assumption that derivatives are small {0{d) = 0{e)). 
From equation 129|) we already know that 0{F) = 
0{d) = 0(e). To avoid obtaining too many terms we 
will limit our analysis here to small velocities, so we also 
postulate 0{U) = 0(e). Note that since we will only ne- 
glect terms of order O(e^) we will not loose any of the 
terms present in the lower order expansion, with the sin- 
gle exception of the V^Q term in equation (|30|l . 

Performing a Taylor expansion to 5th order in the 
derivatives of (|10l) and expressing the fi in terms of 
/o = /o + A, we obtain 

F,; + - tF,) - (r - \)D\f^ - tF,) 
+ {r'~r+l)D^f^^rF,) 

= (40) 

which is the extension of (|25|l to two more orders. Here 
the derivative operator is 

D = dt+v,.V. (41) 

For a stationary profile, advected with velocity U, we 
have the operator identity 

dt + C/.V = 0. (42) 

This simplifies the derivative operator 

D = {v, - U).y. (43) 

We now need to take the zeroth and first order moments 
of this expression to obtain the expressions for the conti- 
nuity and Navier Stokes level equations in the stationary 
advected limit. The expectation is that the continuity 
equation takes the form 

V{pu - pU) = O(e^) (44) 

and the Navier-Stokes level equation will become 

V(P) = 0{e'). (45) 



This then allows us to identify the effective mean fluid 
velocity ii to fifth order and the effective pressure tensor 
P also to fifth order. It is highly desirable, although far 
from obvious, that these quantities be independent of the 
relaxation time r. Otherwise the equilibrium properties 
would be coupled to the transport coefficients. 

This task is cumbersome since it involves up to fifth 
order velocity moments of the equilibrium distribu- 
tion function. We will restrict ourselves here to the 
analysis of a projection of the most common models 
(D2Q7,D2Q9,D3Q15,D3Q19,and D3Q27) to one dimen- 
sion. This projection is the D1Q3 model, a one dimen- 
sional model with three velocities. It is important to note 
that this analysis is entirely sufficient to determine the 
phase-behavior of all the above models. The only issue 
not addressed by this analysis is the isotropy of the mod- 
els. One caveat is that this analysis will only be able 
to make statements of the equilibrium behavior of the 
method, not about the approach to equilibrium. 

A. Analysis of the D1Q3 model 

The D1Q3 model is a one dimensional lattice Boltz- 
mann model with three velocities: Vi — { — 1,0,1}. In 
this model there are only three distinct velocity moments 
at each lattice point, corresponding to the three densities 
fi. The moments of the equilibrium density are 

E./>f =P""+f + A (46) 

Because — Vi, vf — vf etc., all higher order velocity 
moments are given by these first three moments. The 
three moments of the forcing term are given by 

^F, = 0, ^F,w, = F J2f,v^ ^2Fu + -9. (47) 

i i i 

The higher order moments are similarly given by these 
first three moments. 

The algebra involved in calculating the higher mo- 
ments is quite extensive, so I will only outline the method 
here without giving the lengthy intermediate results. 
Summing y^, (|40|) gives an expression for the continuity 
involving all moments of and Fi. Similarly we ob- 
tain such an expression containing all moments by sum- 
ming EiGOl'^i to obtain the momentum equation. We 
can then use the momentum equation to express A in 
terms of the other moments and higher order derivatives 
of A. Then we insert this expression repeatedly into itself 
(similarly to what we did for equation (I24II ) to express 
A completely in terms of the other moments and their 
derivatives. 

This expression for A is then used to eliminate any A 
dependence in the continuity equation. We then use the 
resulting continuity equation to express u in terms of the 
other moments and higher order derivatives of u. This 
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expression for u is then inserted for all but the lowest 
order derivatives of u of the continuity equation. The 
resulting continuity equation has the form 

y{pu -^F- pU) - hjV[FU + \v{pU)] = O(e^) (48) 

which only contains one term of u and no terms with A. 
We use this version of the continuity equation to express 
u in terms of U , p and F and replace all occurrences of u 
in the momentum equation. We then remove all but the 
lowest order occurrences of A in the momentum equation 
and obtain 

F + V [f + A - (r - i)[3™ + V(pC/)] 

+ 4^^^ + T2 - + VI/)] = 0(e5). (49) 

Now we apply these results to our Pressure and Forcing 
methods. 



B. Pressure Method 

For the pressure method we use the moments defined 
in equation and obtain for the continuity equation 

V{pu - pU) - VV[^V(pC/)] = 0(6^). (50) 

Comparing this to H44fl we see that the mean fluid veloc- 
ity is 



1 + i v^- 



(51) 



12 p 



The gradient term is a new correction for the measure- 
ment of the velocity for the pressure method. For the 
momentum equation we obtain 
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FIG. 4: Comparison of the expression for the mean fluid veloc- 
ity from original definition Q and the corrected version H51|l 
for the pressure method in (a). The system was initialized 
with a constant velocity oi U = 0.01. While the correction 
term clearly removes some spurious velocity terms there is 
a more noticeable violation of Galilean invariance which ad- 
vects the gas phase faster than the liquid phase with respect 
to the lattice. In (b) we compare the measured expression 
for the evaporation V[p(?7 — u)] for the original and corrected 
expressions for u and compare it with the theoretical one of 
^Vp for ^ = —0.0002. This simulation was performed with 
= 0.35, e = 1/3 and K = 0.1. 



There are no additional pressure terms for the pressure 
method, which is consistent with the fact that wc found a 
constant input pressure for simulations with the pressure 
method in Figures ^b) andEJb). 

To test the correction for the velocity predicted in H51|) 
we have to consider a situation in which U is not zero. 
We set up a profile that is initiated with a velocity oiU — 
0.01. In Figure ^a) we plot both the standard velocity 
u from ||3 and our new approximation of the true fluid 
velocity u from (|51|) . First we note that the gas phase is 
advected faster than the liquid phase, which is a problem 
of Galilean invariance The gas is advected faster 

than the liquid and a constant evaporation on the leading 
edge of the droplet and a condensation at the trailing 
interface of the liquid leads to a mean interface velocity 
that is less than the imposed mean fluid velocity of 0.01. 
This evaporation and condensation leads to an additional 
velocity ^ of the interface. For such an interface velocity 



we have an additional contribution to the change of the 
density, the rate of evaporation given by 

(9^ -f C/.V)p = -eVp = V[p(C/-u)]. (53) 

So to compare the original dcflnition of ii from equation 
(|27|l with the corrected definition of equation H51I) we plot 
V[p(u — U] for the two definitions of it and ^Vp in Figure 
^b). We see that the corrected version of the velocity 
allows a very good fit with the theoretical expression. 
The best fit for ^ is -0.0002 which corresponds to a 2% 
correction for the observed domain velocity. 

While the Galilean invariance violations are small, in 
this case it is still worth while to correct the Galilean 
invariance violations of the lattice Boltzmann method. 
To do that we need an expansion that does not make the 
assumption u = 0(e). This investigation of higher order 
Galilean invariance violations will be reported elsewhere. 
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C. Forcing Method 

For the forcing method H32|l we obtain the continuity 
equation 



y[pu-^F-U]^0{e^). 



(54) 



This imphes that the mean fluid velocity is given by u = 
pu — F/2 in agreement with the lower order expression 
(|27|l . There are no further corrections to the velocity for 
the forcing method to this level of approximation. 

To evaluate the higher order terms for the pressure in 
equation H49|) we need to insert the specific form of the 
force. We have 



k=0 



y2fc+l 
(2fc+l)! 



(55) 



With this expression for F the momentum equation be- 
comes 



V P-(r 



1 FF 



1. 



-V^F 



(56) 



where we have again expressed some of the higher order 
contributions in terms of the force. We conclude that the 
effective pressure is given by 



peff ^P-{t 



1 FF 
4 P 



1 



V^F + 0{e^ 



(57) 



This explains why the input pressure is not constant for 
the current forcing method. To verify this analytical re- 
sult we can plot both the input pressure P and the pre- 
dicted effective pressure P'^-^'-f for a phase separated sta- 
tionary profile obtained from the forcing version of the 
lattice Boltzmann simulation. The result is shown in Fig- 
ure E] 

Somewhat surprisingly, the pressure P'^f^ is not sim- 
ply a fifth order approximation to the true pressure. In- 
stead, for [/ = 0, this gives a pressure that is constant to 
machine accuracy! To be exact what is constant is the 
discretization 

peff = p _ _ i):^ + iv2(^)p - -v2(^)p. (58) 

^ 4 p 4 12 ^ ^ 

This is, presumably, due to the fact that the higher order 
moments are just repeats of the lower order moments and 
this allows the higher order moments to consist of higher 
order derivatives of the lower order moments. Unfortu- 
nately this exact solution no longer holds ii U 0, be- 
cause the higher order moments contain higher powers of 
U. The fact that we know an exact pressure that we min- 
imize may prove to be very useful for diffusive systems, 
though, for which we have U = 0. 

This may also open the door for an important appli- 
cation of mimetic calculus. If we were able to find some 
discrete gradient operators for which V^P — pV^ p 
holds, we would be able to devise a method that always 
recovers the correct equilibrium behavior. 




30 

lattice position 

FIG. 5: Comparison of the input pressure P from equation 
113711 and the effective pressure P"^^^ of equation 15711 for the 
forcing method for two values of k. We see that the effec- 
tive pressure P'^ff is constant to machine accuracy in both 
cases. This agreement is better than expected since expan- 
sion predicted this pressure only up to fifth order. The critical 
temperature for this simulation is 9c — 0.37 at 6' = 1/3 so that 
true equilibrium pressure is 0.09. 
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FIG. 6: The performance of the corrected Forcing method 
gives results that are identical to the results of the pressure 
method. The interfaces correspond exactly. Therefore the 
pressures and the chemical potentials also agree. Simulation 
parameters are Oc = 0.37, 6 — 1/3 and k = 0.1. 



D. A new forcing method 

Now that we have identified the effective pressure P^ff 
that is constant in the steady state we can identify a 
method that will have a constant input pressure P. Inves- 
tigating equations H48I) and (|49|) we can identify a choice 
of that will allow us to cancel the additional pressure 
terms due to the force. We can amend the original forcing 
method (I22Jl with 



T\E' 



1 FF 
(r- -) — 
^ 4^ p 



+ Y^(vv)V 



(59) 
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FIG. 7: The pressure and the new forcing methods give equiv- 
alent results. We see that the phase-diagram for k = even 
reproduces the non- uniqueness of the solutions. The recovery 
of a universal profile also occurs at the same pace in the new 
forcing and pressure methods, as shown in (b). 

Wc then obtain the continuity equation 

v(^pw-ii^-C/^ =0(6^) (60) 

and the momentum equation 

VP = 0(e5). (61) 

The choice of ^ is not unique regarding the exact choice 
of the correction terms. We could replace (VV)^p with 
— 3V^F and still remain correct up to O(e^). However, 
the above choice again leads to a pressure that is constant 
up to machine accuracy for U — 0, just as we did with 
the pressure method. 

A comparison of the pressure method and the corrected 
forcing method is shown in Figure El The new forc- 
ing method is nearly indistinguishable from the pressure 
method. One small difference between the two methods is 
that an alternating pattern in the pressure does not lead 
to a force, and therefore does not decay. In the simula- 
tion shown in FigurcElthose pressure oscillations have an 
amplitude of about 10^^. However, in some simulations 




FIG. 8: Comparison of the effect of correction terms to the 
forcing method derived by Guo et al. T^] and the current pa- 
per. We use the pressure of equation l|36^ . We observe that 
the agreement with the phase-diagram is improved using the 
correction terms from |15ll . For the new forcing method the 
results would be indeterminate just at the results of the pres- 
sure method in Figure However, including gradient terms 
of ^ with K = 0.1 we recover thermodynamic consistency. 

these oscillations can become large. Also for large veloci- 
ties these oscillations can become unstable and make this 
method less stable than the pressure method at larger ve- 
locities. 

Because the pressure and chemical potential depend 
only on the profile p{x) the new forcing method gives the 
same constant pressure and near constant chemical po- 
tential as the pressure method shown in Figure HJ This 
forcing method is therefore thermodynamically consis- 
tent. It also recovers the analytical phase diagram, as 
shown in Figure |H1 

We also examined the behavior of the forcing method 
in the limit where k ^ 0. This limit may appear tricky 
because the correction terms need to cancel the numerical 
derivatives we uncovered in the expansion. This should 
uncover any higher order terms that we missed in our 
expansion. We find, however, that the method described 
here will even recover the sharp interfaces we observed 
for the pressure method. The non-uniqueness of the so- 
lutions which we found for the pressure method is also 
observed for the pressure method. The recovery of a 
unique solution for larger values of k occurs at the same 
pace as for the pressure method. This is shown in Fig- 
ured This surprising result is due to the fact that the 
pressure P is exactly constant for both methods. 

VI. DISCUSSION 

In this paper we have shown that forcing terms lead to 
non-negligible higher order terms for systems with large 
density gradients. We introduced a new equilibrium anal- 
ysis method, that allowed us to identify correction terms 
up to 5th order. This equilibrium analysis allowed us to 
identified those correction terms. This analysis analysis 
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gave the exact pressure for systems that are not advected 
with respect to the lattice. 

Correction terms similar to, but different from, the cor- 
rection terms we identified have been proposed by Guo 
et al. 01 . Their analysis leads them to conclude that 
for a body force we should choose (in the notation of this 
paper) 

T^=(^r- pFF. (62) 

Comparing this to H59|l we find that this term is different 
in that we predict a factor of 1/4 instead of 1/2 as well as 
an additional double derivative term of p. The correction 
term (|62f) was derived using a multi-scale analysis to sec- 
ond order. We did not understand how these terms could 
be consistently derived with a second order expansion. In 
a Taylor expansion the term ^{pFF) appears only as a 
third order term. It is likely that a second order expan- 
sion would not pick up the VV^p term. If this expansion 
can pick up the pFF term in ^E" consistently, however, 
there cannot be a difference in the pre-factor of this term 
when it is derived by a Taylor expansion method. 

It is therefore important to compare the correction 
term predicted in 15] to our correction term ^ of H59|l . 
We implemented the correction term (|62|l . We then per- 
formed simulations with pressure H36() which does not 
include interfacial terms. We would therefore expect 
a sharp interface. We do, however, observe that this 
method leads to an extended interface, indicating that 
there are additional gradient terms in the pressure not 
accounted for by the Navier-Stokes equation derived by 
Guo et al. 0, equation (18b)]. 

We also measured the phase-diagram for the original 
Forcing method, the corrections proposed by Guo and 
the new forcing method. The results are shown in Fig- 
ure |H1 We see that the correction introduced by Guo et 
al. improves the result for the forcing method somewhat. 
However, there is still no good agreement with the theo- 
retical phase diagram. This is, however, achieved by the 
correction derived in this paper when we include addi- 



tional gradient terms in the pressure tensor H37() . This 
suggests that the prefactor derived in the current work is 
correct. 

There is a longstanding discussion about the suitabil- 
ity of pressure methods for the simulation of non-ideal 
fluids 7, 9]. The criticism relies on the general argument 
that there should be a close correspondence between the 
lattice Boltzmann equation and Kinetic Theory. At this 
level it is a somewhat philosophical argument. And this 
philosophical argument is weakened by the occurrence of 
lattice correction terms that have no analogy in Kinetic 
Theory. A more useful test should be the ability of the 
forcing and pressure methods to recover correct solutions 
in a robust and stable manner. As far as the recovery 
of equilibrium solutions is concerned, we find that there 
is no fundamental difference in the suitability of includ- 
ing non-linear pressure terms into a lattice Boltzmann 
method through a forcing term or a pressure term. 

As a next step in this analysis we need to compare 
the performance of the different methods under Galilean 
transformations. To uncover any term that have not been 
previously discussed |l3 requires us to drop the assump- 
tion of small velocities in our 5th order expansion. Those 
results will be published elsewhere. 

To truly distinguish between the two approaches, how- 
ever, the analysis of dynamic solutions is required. There 
are very few tests in the literature of non-stationary so- 
lution for liquid-gas lattice Boltzmann simulations. But 
such tests will be required to compare the differences be- 
tween pressure and forcing approaches. Examples of such 
simulations include advected fluids undergoing phase- 
separations. For density small compared to the equilib- 
rium densities analytical solutions exist for p{x,t). 

In closing we want to point out that the corrections for 
the forcing method for phase-separated systems do also 
apply for external forces. The same additional pressure 
terms that we identified in (|58|) also occur when a truly 
external force acts on the system. This is particularly 
important to keep in mind when simulating liquid-gas 
systems in the presence of gravity. 
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